home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Eigen-Values and Vectors < prev    next >
Encoding:
Text File  |  1992-09-13  |  3.3 KB  |  107 lines  |  [TEXT/NCII]
  1. ################
  2. # Eigenvalues and Eigenvectors
  3. #
  4. #   program Eigens(theMat,theEigs,theVecs) # Finds theEigs,theVecs from theMat.
  5. #
  6. #
  7. #
  8. #  This text file explains and gives examples of
  9. #  some of the routines in the file All Library Routines, which
  10. #  should be Opened before trying any of these examples.
  11. #################
  12.  
  13. ##############
  14. #  xCOD external routines.  xEIG2 is called by xCOMPLEX_EIG_SPLIT
  15. # and is used for no other purpose. (The code was too big to fit in
  16. # one xCOD, which has a 32k size limit.)
  17. #
  18. xCOMPLEX_EIG_SPLIT
  19.   # xCOD xCOMPLEX_EIG_SPLIT(N:num; matR,matI,eigR,eigI,vecsR,vecsI, w1,w2,w3:RealArray;xCOD xEIG2(dum:num); Err:num); finds eigen values/vectors of a complex matrx. In: real and complex parts matR[1…N,1…N],matI[1…N,1…N]w1[1…n],w2[1…n],w3[1…n] are workspace; xEIG2 does some of the work. Out: eigR,I[1…N, vecsR,I[1…N,1…N]. Err: 0=> success, j=> hit iteration limit on j'th eigen ; an external program.
  20.  
  21. xEIG2
  22.   # xCOD xEIG2(dum:num); only use with xCOMPLEX_EIG_SPLIT.; an external program.
  23.  
  24.  
  25.  
  26. ###############
  27. #  The interface routine:
  28.  
  29. # Eigens calculates the eigenvalues and eigenvectors of
  30. # a square complex (or real) NxN matrix. (N<=1000). 
  31. # The eigens of theMat[1…N, 1…N] are defined by
  32. #      theMat • theVecs[,j] = theEigs[j] * theVecs[,j]
  33. #      for 1<=j<=N.
  34. #
  35.  program Eigens(theMat,theEigs,theVecs) # Finds theEigs,theVecs from theMat.
  36. .   var rm,im,er,ei,vr,vi,i,n,d, err, w1,w2,w3
  37. .  # Input :   
  38. .  #    theMat[1…N, 1…N]  real or complex
  39. .  # Output:
  40. .  #     theEigs[1…N] 
  41. .  #     theVecs[1…N,1…N] complex   
  42. .  #     Err:  0 for no error.
  43. .   begin
  44. .   d=dimension(theMat)
  45. .   if size(d)<>2 then 
  46. .      Print("# •• ERROR : not a matrix.")
  47. .   else if d[1]<>d[2] then 
  48. .      Print("# •• ERROR : not a square matrix.")
  49. .   else
  50. .       begin
  51. .       i = sqrt(-1) # just in case global "i" is different.
  52. .       n = d[1];
  53. .       w1=1…n;  w2=w1; w3=w1;
  54. .       rm = real(theMat)
  55. .       im = imaginary(theMat)
  56. .       er=1…n; ei = er; vr=rm;vi=rm; err=0;
  57. .       xCOMPLEX_EIG_SPLIT(n,rm,im,er,ei,vr,vi, w1,w2,w3, xEIG2,err)
  58. .        if err<>0 then Print("# •• ERROR in Eigens = "+err)
  59. .        theEigs = er + i*ei
  60. .        theVecs = vr+i*vi
  61. .        end
  62. .    end
  63. .    
  64.  
  65. ########################
  66. # For example, 
  67.  
  68. a = [0,0, 1; 0,1,0; 1,0,0]    # Define a matrix "a".
  69.   a[1…3,1…3] = [[0, 0, 1], 
  70. .                         [0, 1, 0], 
  71. .                         [1, 0, 0]]  .  
  72.  
  73. Eigens(a,aVals,aVecs)           # Find it's eigens.
  74.  
  75. avals    # The three eigenvalues:
  76.   aVals[1…3] = [1, -1, 1] 
  77.  
  78. avecs   # and the eigenvectors:
  79.   aVecs[1…3,1…3] = [[0.70711, -0.70711, 0], 
  80. .                                [0, 0, 1], 
  81. .                                [0.70711, 0.70711, 0]]  .  
  82.  
  83. # The egienvalues act like the entire matrix "a" on the columns
  84. # of the eigenvector matrix. For example, with the 2nd column:
  85.  
  86. a • aVecs[,2]                           # This is the same as
  87.   [0.70711, 0, -0.70711] 
  88.  
  89. aVals[2] * aVecs[,2]               # this.
  90.   [0.70711, 0, -0.70711] 
  91.  
  92.  
  93. # Note that each eigenvector is normalized so that it has
  94. # length 1.  (Since the vectors may be complex, ust † (opt-t) to conjugate
  95. # the second.)
  96. aVecs[,2] • aVecs[,2]†
  97.   1 
  98.  
  99.  
  100. ##########
  101.  
  102. # Remember,  "•" (opt-8 on US keyboard) is matrix "dot" multiplication,
  103. #                    "*" is "normal" element-by-element multiplication.
  104. #  (See Multiply Options in the Other menu.)
  105.  
  106.  
  107.